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Abstract 

A  suite  of  astrodynamic  software  tools  has  been  developed  for  associating  un-catalogued  analyst  satellite  debris  with 
historical  launches,  breakups,  and  anomalous  events.  A  semi-analytical  orbit  integrator  is  at  the  heart  of  the  tool 
suite.  To  associate  currently  tracked  analyst  debris  with  events  occurring  over  30  to  45  years  ago  requires  an  orbit 
integrator  that  predicts  orbit  plane  ascending  node  values  accurately,  with  errors  no  more  than  a  few  degrees  over 
thousands  of  degrees  of  node  precession.  Such  an  integrator,  SGPE  (Special  General  Perturbations  Ephemeris),  has 
been  developed  and  proven  to  be  extremely  accurate.  Orbit  perturbations  that  are  included  in  the  tool  are  daily 
atmospheric  density  corrections,  atmospheric  rotation  variations,  high-order  lunisolar  gravitational  effects,  earth 
tidal  effects,  solar  radiation  pressure  perturbations,  and  high-order  geopotential  zonal  contributions.  All  these 
perturbations  are  computed  with  analytical  equations  evaluated  on  a  half-day  increment  to  achieve  a  very  fast  semi- 
analytical  integrator.  Also  part  of  the  tool  suite  is  a  non-linear  least  squares  software  package  used  to  determine 
average  ballistic  coefficient  values  obtained  from  using  years  of  mean  element  semi-major  axis  values.  Once  the  B 
value  has  been  fitted  the  SGPE  integrator  is  used  to  integrate  backwards  from  epoch  to  1958,  generating  an 
ephemeris  of  mean  element  sets  every  5  days.  These  elements  are  then  compared  with  mean  elements  from  every 
launch  and  breakup  since  the  1960s.  Probabilities  of  possible  associations  are  computed  from  tables  based  on  over 
550  validation  cases  using  historical  mean  elements  from  catalogued  debris.  Associations  of  hundreds  of  analyst 
satellite  debris  with  historical  events  have  already  been  identified,  and  cataloging  of  the  analyst  satellites  is  currently 
underway.  Examples  are  shown  of  breakup  associations,  and  of  associating  currently  tracked  debris  with  previously 
lost  satellites. 

1.  Introduction 

Numerous  attempts  have  been  made  to  associate  currently  tracked  un-cataloged  analyst  debris  with  historical 
launches  and  breakups.  Prior  to  1980  there  have  been  numerous  breakups,  as  shown  in  Fig.  1.  The  backward 
integration  must  span  at  least  25  to  45  years  to  associate  with  very  early  events.  To  associate  debris  with  a  particular 
launch  or  breakup  the  orbit  plane  orientation  needs  to  be  accurately  correlated.  This  means  that  the  orbit  inclination 
and  right  ascension  of  the  ascending  node  need  to  be  accurately  computed.  The  problem  is  to  define  the  right 
ascension  of  the  ascending  node  to  an  accuracy  of  only  a  few  degrees  after  the  very  long  integration  span.  Table  1 
shows  the  amount  of  precession  of  the  node  over  a  40  year  span  as  a  function  of  the  orbit  inclination.  From  12,000 
to  over  65,000  degrees  of  precession  must  be  accurately  computed  for  the  great  majority  of  the  events.  In  previous 
association  attempts  the  integration  has  been  done  using  SGP4  or  SP,  the  special  perturbation  integrator.  Using  the 
analytical  mean  element  SGP4  propagator  means  fast  long  term  integrations  are  possible,  but  the  predictions  of  the 
node  are  very  inaccurate  because  of  a  lack  of  in-depth  perturbation  computations.  Fig.  2  shows  examples  of  four 
different  long  term  SGP4  propagations  for  orbits  of  different  inclinations  and  altitudes.  Errors  greater  than  1000 
degrees  can  occur  using  this  method.  Using  a  special  perturbation  integrator  increases  the  accuracy  of  the 
integration  but  drastically  increases  the  computer  time  for  the  integration.  However,  as  will  be  shown  later,  using 
special  perturbations  still  does  not  accurately  account  for  all  the  necessary  perturbations  to  achieve  the  desired 
accuracy  of  the  node  predictions. 
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Breakup  Events  by  Year  for  Current  Pieces  in  Orbit  >=  5 


Fig.  1.  Breakup  events  by  year  for  breakups  with  more  than 
5  pieces  currently  in  orbit. 


Inclination 

RA  Precession 

Degs  after  40  Years 

(deg/day) 

50° 

-4.60 

-67200° 

70° 

-2.30 

-33600° 

82° 

-0.82 

-12000° 

90° 

0.01 

+  - 

98° 

0.95 

13880° 

Table  1.  Precession  of  the  right  ascension  (RA)  of  the  orbital  node  as  a  function 
of  orbit  inclination,  with  the  number  of  degrees  of  precession  after  40  years  listed. 


A  RA  Values  Using  J2,  J3,  Ndot  Prediction 


Fig.  2.  Error  in  right  ascension  of  ascending  node  (ARA)  for 
selected  satellites  using  SGP4  orbit  perturbations. 


2.  Semi-analytical  Orbit  Integrator  SGPE 


To  achieve  an  accuracy  of  only  a  few  degrees  in  the  precession  of  the  orbit  node  requires  in-depth  perturbation 
computations  to  cover  40  years  or  more  of  time.  A  one  step  propagator  like  SGP4  would  require  exhaustive 
analytical  equations  accounting  for  all  the  variations  of  the  orbit  plane  as  a  function  of  time.  Using  special 
perturbation  techniques  with  detailed  small  step  integration  would  require  many  hours  of  computer  time.  The 
solution  for  this  problem  is  to  use  a  semi-analytical  integration  technique  with  a  step  size  on  the  order  of  a  day  that 
accurately  accounts  for  all  the  perturbations.  The  semi-analytical  integrator  SGPE  was  originally  developed  in  the 
1970s  to  predict  long  term  orbit  perturbations  of  high  eccentricity  satellites.  Only  first-order  perturbations  of 
geopotential  resonance  and  lunisolar  gravitational  effects  were  considered  for  the  initial  model.  It  was  determined 
that  a  half-day  step  size  would  be  sufficient  to  achieve  accurate  integrations.  Recently  SGPE  was  extended  to 
accurately  predict  node  precession  by  adding  high  order  drag  computations,  increased  lunisolar  gravitation  effects, 
earth  tides  changes,  solar  radiation  pressure,  and  increased  geopotential  zonal  perturbations. 

From  previous  drag  analysis  (Bowman  [1])  it  was  determined  that  using  any  current  empirical  atmospheric  model  in 
any  integrator  would  result  in  inaccurate  modeling  of  the  atmosphere,  resulting  in  inaccurate  atmospheric  drag 
computations.  However,  adding  a  daily  density  correction  to  the  model  would  greatly  reduce  the  errors  in  the  drag 
computations.  Fig.  3  shows  the  results  of  the  ballistic  coefficient  solutions  based  on  using  the  Jacchia  70  [2]  model 
with  and  without  a  daily  density  correction  based  on  observed  density  variations.  Special  perturbation  orbit  fits 
were  used  for  the  B  solutions.  Bowman  [1]  describes  the  method  of  using  calibration  satellites  at  different  altitudes 
to  compute  daily  temperature  (i.e.  density)  corrections.  This  method  was  used  to  compute  daily  density  corrections 
as  a  function  of  altitude  from  1962  through  2009.  The  satellite  00060,  used  for  the  B  computations  in  Fig.  3,  is  a 
blunt  double-cone  shape  which  behaves  like  a  sphere  for  drag  computations.  Thus,  the  errors  in  the  B  values  are  a 
direct  reflection  of  the  errors  in  the  atmospheric  density.  Fig.  3  shows  that  the  density  errors  are  significantly 
reduced  from  over  17%  one  standard  deviation  to  just  under  6%  by  applying  an  observed  daily  temperature  (density) 
correction  to  the  model  during  the  B  orbit  fits. 


Satellite  00060,  Explorer  8,  DB/BTrue  Values  at  400  km 


Fig.  3.  Ballistic  coefficient  B  error  (%)  from  orbit  differential  corrections  as  a  function  of  year 
(1970  to  2000)  using  the  Jacchia  70  model  with  and  without  daily  temperature  corrections  (dTc). 


King-Hele  [3]  analytical  drag  computations  at  the  orbit  perigee  location  were  selected  for  use  in  SGPE,  coupled  with 
the  CIRA  72  [4]  atmospheric  density  model  using  the  daily  temperature  corrections.  The  atmospheric  drag  on  the 
orbit  is  recomputed  every  SGPE  half-day  integration  step.  The  CIRA  72  model  was  selected  because  at  every 
integration  step  the  density  and  atmospheric  molecular  constituent  amounts  are  computed  from  an  altitude 
integration  of  the  diffusion  equations.  This  is  necessary  because  the  drag  coefficient  is  a  function  of  molecular 
number  densities  (  Bowman  [5]  )  which  varies  with  altitude  and  solar  conditions,  as  displayed  in  Fig.  4.  Therefore, 
the  drag  coefficient  value  needs  to  be  reevaluated  at  every  SGPE  integration  step  for  the  drag  computations.  Also, 
the  atmospheric  rotation  rate  needs  to  be  used  to  account  for  the  drag  effect  on  the  orbit  plane  orientation.  King- 
Hele  [6]  determined  that  the  zonal  winds  are  variable  as  a  function  of  altitude,  so  these  effects  were  also  added  into 
the  SGPE  drag  computations.  For  low  eccentricity  orbits  (e  <  0.05)  it  was  determined  that  the  day-night  diurnal 
density  variations  were  too  large  to  use  the  King-Hele  equations  evaluated  just  at  the  perigee  point  of  the  orbit. 
Therefore,  additional  development  in  SGPE  was  undertaken  to  account  for  the  diurnal  effects  for  low  eccentricity 
orbits. 
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Fig.  4.  Change  in  ballistic  coefficient  B  values  as  a  function  of  year  from  changes  in  the 
drag  coefficient  Cd  due  to  changes  in  atmospheric  constituent  densities  at  high  altitudes 


Fig.  5  shows  the  results  of  the  backward  SGPE  integration  from  2010  to  1968  for  satellite  02826,  a  calibration 
sphere,  with  and  without  applying  a  daily  temperature  (density)  correction.  All  the  other  orbit  perturbations, 
discussed  below,  were  included  in  the  42  year  integrations.  The  SGPE  integration  with  the  daily  corrections 
compares  almost  precisely  with  the  historical  operational  mean  elements,  while  the  SGPE  integration  without  the 
daily  corrections  is  in  error  in  orbital  period  by  about  0.25  minutes  after  the  42  year  integration.  Even  though  this 
period  difference  appears  small  it  is  enough  to  cause  a  large  error  in  the  precession  of  the  node  as  displayed  in  Fig. 
6.  The  error  in  the  SGPE  integration  with  the  daily  corrections  is  just  less  than  6  degrees  after  42  years,  while  the 
integration  without  the  daily  corrections  is  in  error  by  almost  80  degrees  after  the  42  year  backward  integration. 
Thus,  using  the  daily  density  corrections  is  absolutely  vital  to  obtaining  accurate  node  precession  predictions  over 
long  periods  of  time. 


02826  Period  Values  Q  =  -850  km  e  =  0.001  i  =  69.9  B  =  0.21 


Fig.  5.  Variation  in  orbital  period  of  sphere  02826  as  a  function  of  backward  integration  in  days 
since  1950,  with  and  without  daily  temperature  corrections,  plotted  with  operational  mean  elements. 
Q  is  the  perigee  height  for  this  near  circular  orbit. 


02826  dRA  Values  Q  = -850  km  e  =  0.001  i  =  69.9  B  =  0.19 


Fig.  6.  Error  in  right  ascension  of  node  of  sphere  02826  as  a  function  of  backward 
integration  in  days  since  1950,  with  and  without  daily  temperature  corrections  (dTC). 


For  lunisolar  gravitational  perturbations  the  original  SGPE  used  first-order  theory  equations  obtained  from 
integration  of  the  first  term  of  the  disturbing  function.  These  analytical  equations  accounted  for  secular  and  long- 
periodic  perturbations  of  the  orbital  elements.  These  perturbations  were  then  extended  in  the  updated  SGPE  from 
Kozai  [7]  who  used  the  first  5  terms  of  the  disturbing  function  to  compute  analytic  equations  of  the  orbital  element 
variations  as  a  function  of  the  polar  coordinates  of  the  sun  and  moon.  These  equations  are  ideal  for  use  in  a  semi- 
analytical  integrator  with  a  half-day  step  size,  where  the  perturbational  effects  are  continually  updated  each  step  with 
new  sun  and  moon  positions.  Kozai  [7]  also  included  the  perturbations  due  to  the  tidal  deformation  of  the  earth, 
expressed  again  as  first-order  analytical  equations.  Fig.  7  displays  the  SGPE  integration  of  satellite  00900,  a 


calibration  sphere  at  1025  km  altitude.  The  integration  is  with  and  without  earth  tides,  and  is  compared  to  the 
historical  operational  mean  elements.  The  precession  of  the  node  is  highly  dependant  on  the  orbital  inclination,  and 
errors  in  the  inclination  result  in  large  errors  in  the  precession  of  the  node.  The  SGPE  integrations  in  Fig.  7  include 
all  other  perturbations  discussed  in  this  paper.  At  an  inclination  of  precisely  90  degrees  there  should  not  be  any 
node  precession  cause  by  geopotential  zonal  effects  (see  below).  However,  the  solid  earth  tidal  effects  move  the 
orbit  plane  away  from  the  precise  90  degree  inclination  orbit,  and  thus  cause  the  node  to  start  precessing  due  to 
zonal  perturbations.  Fig.  8  shows  the  results  of  the  nodal  precession  mainly  due  to  these  tides.  Accounting  for  the 
tides  in  SGPE  produces  a  45  year  ephemeris  that  almost  precisely  mirrors  the  historical  observed  variations  in  the 
inclination.  However,  very  large  errors  occur  in  the  integration  by  not  including  the  earth  tides  effects. 
Additionally,  general  precession  and  nutation  effects  were  included  in  SGPE  ,  both  of  which  are  caused  by  the 
gravitational  attraction  of  the  sun  and  moon  on  the  equatorial  bulge  of  the  earth  (Cooley  [8]  and  Newton[9] ). 


00900  i  Values  Q  =  1025  km  e  =  0.002  i  =  90.2  B  =  0.24 


Fig.  7.  Variation  in  inclination  of  sphere  00900  as  a  function  of  backward  integration  in 
days  since  1950  with  and  without  tidal  effects,  plotted  with  operational  mean  elements. 


00900  RA  Values  Q=  1025  km  e  =  0.002  i  =  90.2  B  =  0.24 


Fig.  8.  Right  ascension  of  node  of  00900  as  a  function  of  backward  integration  in  days 
since  1950,  with  and  without  tidal  effects,  plotted  with  operational  mean  elements. 


The  geopotential  zonal  perturbations  produce  the  main  effect  of  the  precession  of  the  orbital  plan.  The  very  large 
node  changes  listed  in  Table  1  are  the  result  of  this  perturbational  effect.  These  variations  are  represented  by  the 
even  zonal  coefficients  of  the  earth’s  geopotential  field.  SGP4  uses  only  the  first  two  major  coefficients  J2  and  J3  for 
the  secular  and  long-periodic  variations  in  the  orbital  elements.  However,  for  an  accurate  long-term  orbit  integration 
many  more  higher  order  terms  are  needed.  Mueller  [10]  computed  analytical  equations  for  these  perturbations  with 
all  inclusive  terms  through  Ji0,  including  perturbations  on  the  node,  inclination,  eccentricity,  and  argument  of 
perigee.  These  equations  were  included  in  SGPE.  The  equations  were  further  extended  to  include  J!2  and  .1 1 4  from 
Kozai  [11]  for  additional  accuracy.  Finally  the  J2  squared  term  from  King-Hele  [12]  was  included  in  SGPE  for  the 
node  perturbation  to  capture  the  use  of  Kozai  type  mean  elements  used  in  the  SGPE  integration. 

The  last  SGPE  included  perturbation  of  the  node  affecting  low  earth  orbits  is  solar  radiation  pressure.  Koskela  [13] 
developed  a  complete  theory  for  handling  the  effects  of  solar  radiation  pressure  on  satellite  orbits.  The  effect  of  the 
satellite  being  eclipsed  by  the  Earth  was  considered,  as  well  as  the  case  when  the  satellite  is  continuously  exposed  to 
sunlight  for  many  orbits.  The  eccentric  anomaly  was  used  as  the  independent  variable  in  the  analytical  approach, 
with  perturbational  equations  being  obtained  in  closed  form.  The  key  to  obtaining  an  accurate  solar  radiation 
perturbation  is  determining  an  accurate  radiation  pressure  coefficient.  The  non-linear  least  squares  solution  for  B 
below  discusses  the  solar  radiation  pressure  coefficient  determination. 

3.  Non-Linear  Least  Squares  Solution  for  B  Value 

A  major  key  to  accurate  long  term  predictions  is  an  accurate  atmospheric  drag  computation,  which  means  using  an 
accurate  ballistic  coefficient  B  value  with  accurate  density  values.  A  non-linear  least  squares  model  was  developed 
for  use  with  the  SGPE  integrator.  This  model  uses  multiple  years  of  mean  element  semi-major  axis  values  to  solve 
for  an  average  B  value  and  an  initial  mean  motion  value.  The  great  majority  of  satellite  debris  has  B  variations  due 
to  frontal  area  variations  resulting  from  the  non-spherical  shapes  of  the  pieces.  Therefore,  it  is  necessary  to  obtain 
the  best  average  value  possible  for  long  term  predictions.  Obtaining  an  accurate  B  average  value  assumes  that  the 
rotation  and  precession  of  the  rotation  axis  of  the  pieces  will  have  periodicities  shorter  than  the  fit  span  used  for  the 
average  B  determination.  More  than  550  known  catalogued  pieces  of  debris  were  used  in  the  validation  process  to 
determine  integration  accuracy  as  a  function  of  the  B  fit  span  (in  years),  the  magnitude  of  the  B  value,  the  altitude  of 
the  orbit,  and  the  length  of  the  SGPE  integration.  One  sigma  standard  deviation  error  tables  were  computed  from 
the  validation  cases,  and  these  are  used,  as  described  below,  to  determine  boundary  values  for  assigning  event 
association  probabilities. 

The  non-linear  least  squares  program  is  designed  to  process  at  least  3  iterations  until  the  B  solution  reaches  a 
consistent  value.  Numerical  partials  are  computed  for  iterating  on  the  initial  mean  motion  and  B  values.  Following 
the  B  solution  the  left  over  residuals  in  inclination,  right  ascension,  and  argument  of  perigee  are  then  least  squares  fit 
to  account  for  any  remaining  unmodeled  perturbation  effects.  These  new  terms  are  then  applied  as  additional 
perturbations  in  the  SGPE  backward  integration  to  1958. 

Once  the  long-term  average  B  value  has  been  determine  the  SGPE  integrator  is  used  to  generate  2 -line  mean 
elements  sets  every  1  to  5  days  back  to  1958.  This  mean  element  “ephemeris”  is  then  compared  with  the  mean 
elements  from  initial  launch  and  breakup  elements  of  the  parent  satellites  and  rocket  bodies.  Databases  of  mean 
elements  were  built  using  the  first  30  days  of  elements  from  every  launch,  and  elements  from  20  days  before  and 
after  every  breakup  that  occurred  since  the  1960s.  All  these  launch  and  breakup  element  sets  (over  200,000  for 
launches  and  over  230,000  for  breakups)  are  then  compared  to  the  element  sets  generated  in  the  SGPE  ephemeris 
file.  The  probability  of  association  of  a  SGPE  ephemeris  element  set  with  a  launch  or  breakup  element  set  is  based 
upon  the  integration  error  obtained  from  tables  of  validation  data  described  above.  This  association  computer  run 
(ELMCOMP)  only  takes  a  few  seconds  of  computer  time  to  complete.  The  longest  computer  time  required  for  this 
entire  process  is  for  the  SGPE  integration,  which  takes  approximately  1  second  per  integration  year. 

4.  Analyst  Debris  Association  Examples 

Following  an  association,  the  launch  or  breakup  parent  element  sets  and  the  SGPE  ephemeris  element  sets  are 
plotted.  Fig.  9  is  an  example  of  associations  of  several  analyst  debris  with  2  separate  breakups  occurring  in  late 
1992  and  early  1993.  To  remove  the  very  large  nodal  precession  that  occurred  over  the  last  two  decades  the  element 


set  node  values  of  the  parent  satellite  22285  are  used  as  a  reference.  The  node  differences  from  the  parent  values  of 
the  analyst  debris  are  then  computed  and  this  “relative”  right  ascension  of  the  ascending  node  is  plotted  for  each 
debris  piece.  Fig.  9  shows  very  high  correlations  (intersection  of  the  relative  nodes)  of  analyst  debris  81535  and 
89173  with  the  breakup  of  satellite  22285  in  late  1992.  Also,  analyst  debris  81215  and  89115  shows  very  high 
correlations  with  the  breakup  of  22566  occurring  in  early  1993.  All  of  the  elements  of  the  analyst  debris  were 
generated  from  an  SGPE  integration  going  backward  from  the  mid  2000s.  This  method  of  first  computer  association 
and  then  verification  by  plotting  have  been  used  to  identify  many  hundreds  of  pieces  associated  with  historical 
launches  and  breakups. 


Fig.  9.  Associations  of  analyst  satellite  nodes  with  2  separate  breakups  occurring  in  1992  and 
1993.  The  relative  right  ascension  of  the  ascending  nodes  are  plotted  as  a  function  of  date. 


In  addition  to  associations  with  historical  events  this  association  method  can  be  used  to  recover  lost  satellites.  Fig. 
10  displays  an  example  of  identifying  a  newly  tracked  piece  of  debris  with  a  previously  lost  satellite.  Satellite 
23293  is  a  moderately  eccentric  small  piece  of  SL-12  debris  with  an  apogee  of  -19,000  km  and  a  perigee  height 
varying  between  500  km  and  900  km.  The  argument  of  perigee  of  the  orbit  has  a  very  slow  precession,  and  this  fact 
coupled  with  the  varying  perigee  height  produced  conditions  such  that  the  debris  was  not  able  to  be  tracked  after 
mid  2005.  The  green  curve  in  the  figure  shows  the  tracking  from  breakup  to  2005.  The  non-linear  least  squares  B 
solution  was  obtained  with  the  elements  prior  to  loss  of  tracking,  and  the  SGPE  integrator  was  used  to  generate  an 
ephemeris  of  element  sets  from  2005  through  2010  (red  curve).  This  element  set  ephemeris  was  then  compared  with 
known  analyst  debris  and  a  high  correlation  was  obtained  with  debris  piece  83984  that  was  only  tracked  since  the 
beginning  of  2010.  Fig.  10  displays  the  correlation  in  inclination,  node,  apogee  and  perigee  heights,  period,  and 
eccentricity,  and  leaves  no  doubt  that  83984  is  indeed  the  lost  23293  piece.  This  example  illustrates  the  use  of 
SGPE  not  just  for  historical  event  associations  but  also  for  the  recovery  of  lost  satellites. 
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Fig.  10.  Mean  element  plots  of  satellite  23293  (green)  up  until  time  of  loss  (2005), 

SGPE  predicted  elements  (red),  and  analyst  satellite  83904  elements  (blue). 

5.  Conclusions 

It  has  been  demonstrated  that  currently  tracked  un-catalogued  analyst  satellites  can  be  associated  with  historical 
launch,  breakup,  or  other  anomalous  events  that  have  occurred  as  far  back  as  the  1960s.  A  fast  semi-analytical 
integrator  SGPE  has  been  developed  using  analytical  equations  to  account  for  the  great  majority  of  perturbations  of 
the  orbit  plane.  Additionally  a  fast  accurate  method  of  determining  the  average  ballistic  coefficient  has  been 
developed  using  a  non-linear  least  squares  program  fitting  standard  2-line  semi-major  axis  mean  elements. 
Associations  of  currently  tracked  satellite  debris  with  historically  lost  satellites  are  now  also  possible.  Cataloging 
the  currently  thousands  of  pieces  of  analyst  debris  with  accurate  international  designators  is  now  possible. 
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